GENERAL CIRCULATION EXPERIMENTS WITH 
A TWO-LEVEL QUASI -GEOSTROPHIC MODEL 
INCLUDING THE NON-LINEAR INTERACTION 
BETWEEN A SINGLE WAVE IN THE ZONAL 
DIRECTION AND THE MEAN FLOW 



Frank H. Taylor 



United States 
Naval Postgraduate School 




THESIS 



GENERAL CIRCULATION EXPERIMENTS 
WITH A TWO-LEVEL QUASI-GEOSTROPHIC MODEL 
INCLUDING THE NON-LINEAR INTERACTION BETWEEN 
A SINGLE WAVE IN THE ZONAL DIRECTION AND THE MEAN FLOW 



by 



Frank H. Taylor 



April 1970 



Thu document kcu> been approved ^on. pubtic Ke- 
Zecu>e and 6 ate; iti> dLuVubation. u unlimited. 



J133729 




General Circulation Experiments 
With a Two-Level Quasi-Geostrophic Model 
Including the Non-Linear Interaction Between 
A Single Wave in the Zonal Direction and the Mean Flow 



by 



Frank H, ^Taylor 

Lieutenant Commander, United States Navy 
B.S f E E,, Purdue University, 1961 



Submitted in partial fulfillment of the 
requirements for the degree of 



MASTER OF SCIENCE IN METEOROLOGY 



from the 



NAVAL POSTGRADUATE SCHOOL 
April 1970 



ABSTRACT 



A long-period forecast is made utilizing a two-level quasi- 
geostrophic model. The model includes friction and heating which is 
a linear function of y. The model was further simplified by restrict- 
ing the disturbance to one wave in the zonal direction. Experiments 
were performed with two distances between the walls. In the case of 
the longer separation, a solution with appreciable time fluctuations 
was obtained with the largest fluctuation being investigated in more 
detail. The smaller separation revealed a traveling wave of constant 
amplitude similar to that observed in certain dishpan experiments. 
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I. INTRODUCTION 



Knowledge of the general circulation has been made possible with 
the advent of high-speed large-capacity electronic digital computers. 

In a classic experiment, Phillips [1956] performed a long-period 
numerical forecast with a two-level quas i-geos trophic model in a 
zonally periodic channel on a (3 plane. Friction and heating were 
included in his equations with the latter as a linear function of 
latitude. He introduced a finite amplitude initial disturbance into 
a baroc linica 1 ly unstable zonal current. Major deficiencies were the 
error from the |3 plane and the geostrophic approximations, the latter 
making it impossible to account for non-geos trophic dynamics . Further- 
more , due to computational errors, he did not achieve a statistically 
steady state. However, Phillips did initially demonstrate the feasi- 
bility of numerical simulation of the atmosphere. 

General circulation experiments have progressed in sophistication 
and complexity since Phillips’. Smagorinsky [1963] formulated a two- 
level baroclinic primitive equation model by using the Eulerian 
equations of motion within a spherically zonal strip. Continued 
sophistication of the experiments resulted in a nine-level model 
utilizing the primitive equations of motion [Smagorinsky, Manabe 
and Holloway, 1965]. This study included a resolution of boundary 
layer fluxes and radiative transfer involving ozone, carbon dioxide 
and water vapor. A continuation of this study was performed for the 
tropics which was successful in simulating the tropical convergence 
zone and the thermal structure [Manabe and Smagorinsky, 1967]. Mintz 
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[1964] conducted an experiment in which the primitive equations of motion 
were developed for a two-level model of the atmosphere with the heating 
and friction terms retained. This model further brought in the effects 
of land-sea contrasts in heating. Nitta [1967] utilized a twenty-level 
model for computing vertical distribution of the geopotential flux and 
conversion of eddy potential energy and eddy kinetic energy due to non- 
stationary disturbances. A model which has the distinguishing feature 
of using height rather than pressure was developed [Kasahara and 
Washington, 1967] and hydrostatic equilibrium was maintained in the 
system. Some advantages of this model are the prognostic equations 
have a simpler form than those of the a system and the lower boundary 
conditions may be easily formulated. The integration of the nine-level 
primitive equation model was further extended by the box method [Kurihara 
and Holloway, 1967] which conserved the conservat iona 1 properties of the 
original equations. Recently, the effect of hydrology of the earth's 
surface was incorporated in the nine-level primitive equation model 
[Manabe , 1969]. The scheme involved the prediction of water vapor in 
the atmosphere and the prediction of soil moisture and snow cover. The 
numerical integrations were performed for the annual mean distribution 
of solar insolation. 

The increased sophistication and complexity of recent models have 
also produced complex results which make interpretation difficult. 
Therefore, in this study, the complexities have been kept to a minimum 
while still maintaining a description of the basic dynamic processes. 

When the effects of heating and friction are included, the two-level 
quasi-geostrophic model is the simplest model capable of describing the 
general circulation. 
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The model is further simplified in that the disturbance is restricted 
to a single wave in the zonal direction. Phillips* [1956] experiments 
and Smagorinsky * s [1963] experiment which contained no forcing function 
in the east-west direction, displayed a single predominate wave number 
which contained much of the disturbance energy. The y-structure of the 
disturbance and the mean flow are calculated accurately with a sufficient 
number of grid points. Some general circulation experiments have been 
made with even simpler models which represent the y-variation of all 
quantities with a few terms of a Fourier series [Bryan, 1959; Lorenz, 
1962, 1963; Young, 1966], These models do not provide an adequate 
description of the barotropic interaction between the disturbance and 
the mean flow. Consequently, this model appears to be the simplest 
model which can describe the main mechanisms of the general circulation 
of the atmosphere. 
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II. THE FORECAST EQUATIONS 



A simple two-level model is constructed by dividing the entire 
atmosphere into four layers of constant pressure differential, ^ 
(fig. 1), numbered 0 to 4 from top to bottom. Assume vertical motion 
to be zero at the top of the atmosphere, while 
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Fig. 1. Two- level model used for prediction 



the vertical motion term at the earth's surface is approximated by 

CO ~ -pgw 4> (2.1) 

Charney and Eliassen [1949] used the Ekman theory to derive an 
expression for w^ which is 

2 A 1/2 

“ 4 ” J ( — sin 2 a C 4 . (2.2) 

o 

where A^ is the vertical eddy viscosity and a, the inflow angle. The 
surface geostrophic vorticity is approximated by 

C 4 = C 3 = V 2 Y 3 - (2.3) 

Following the development of Thompson [1961], but with different 
notation, begin with the quas i-geos trophic vorticity equation 

V 2 r + k x VfV(V 2 Y) + e 0 g - f 0 | = A h vV (2.4) 
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Apply this equation at levels 1 and 3 giving 



dY, 



If v \ + lk x v V v(v2t i> + B o ST - £ o ff - a h v \' (2 - 5 > 



S V 2 Y 3 + ik x VY 3 -V(V 2 Y 3 ) + 0, 



dY, (<u, -cd ) 

+ f 



a h v Y 3 . 



( 2 . 6 ) 



dt v t 3 ^ a vi 3' vvv t 3 ^ - p q ^ — Q Ap 

If it is desired to stop at some intermediate point, say the tropopause, 
smaller layers could be uti lized rather than to include the entire atmosphere. 

Next consider the quasi-geostrophic first law of thermodynamics in 
the form 



+ Ik * VY-V (I*) + f- a, - -^0. 



f o P 



where 



a = 



R t r St 

s 



2 

P § 



__s + 

<3z C 



]• 



(2.7) 



(2.7a) 



and T g is the temperature from the standard atmosphere. When this 
equation is applied at level 2, the result is 

Ap cm 



N (Yi + Yo) “F 

(y^) + |k x v — r*~ ' V <vv — t 



H ^ 

T ^2 



( 2 . 8 ) 



Define the following quantities 



m 



Y = 

Y T 



Y i + Y 3 



Y - Y 

T 1 t 3 



(2.9) 

(2.10) 



which implies 



Yi = Y m + Y t , 
Y 3 = Y m - Y t . 



( 2 . 11 ) 

( 2 . 12 ) 



Here is proportional to the layer thickness and is therefore a 
measure of the mean temperature. Using these definitions, add (2.5) 
and (2.6) and divide the result by 2, obtaining 
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(2.13) 





For the second forecast equation, subtract (2.6) from (2.5) and 



eliminate usi^ng (2.8) 





(2.14) 



where 



2 




(2.14a) 



These are the prediction equations for the model. At this point 
we restrict ourselves to a disturbance of wave number k and the mean 
flow. The fields may be defined as follows 



where k is the x wave number, A through D are Fourier amplitudes of the 

disturbance and E and F are the zonal mean fields. Now substitute 

expressions for Y and Y m into (2.13), separate the various terms, 
m T 

neglecting all terms with wave number 2k or higher. Equating coeffici- 
ents of the cosine terms gives 



Y = E(y,t) + A(y,t) cos kx + B(y,t) sin kx, (2.15) 
m 

Y t = F(y,t) + C(y,t) cos kx + D(y,t) sin kx, (2.16) 





B + ~ 



§ D )] - e o Bk - K - K 2 ) (A-C) + 
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For the sine terms the result is 



S fb^B rn 2^ r SE 5 A SF d C , S E 
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oy oy 

Repeating the procedures for (2.14), the equations are 



(2.18) 
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and 
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(2.20) 
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Equating coefficients of terms independent of x in (2.13) gives 



a a 2 E _ k a r a 2 B 

dt ay 2 - 2 a y L ^ y 2 



B^f +C 

5y 



d 2 D 

Sy 2 



»f§]- 

dy 



K 



c>_ fh £ . . * iL E " N i 

dy V<3y byJ n <3y V^ y 3y‘ 



(2.21) 



Similarly, coefficients indepdendent of x in (2.14) results in 

rb 2 kb r, a 2 D _ a 2 c n2 - 

vtt “ ^ J ^ " 2 a; L A 72 - B tt + c 



% 



dy 



dy 



£2 . D £± + 
8y 2 V 



where 




( 2 . 22 ) 



(2.23) 



The prediction equations are now in the form to be used in the model. 



14 



III. FINITE DIFFERENCE SCHEME AND BOUNDARY CONDITIONS 



The finite difference scheme used is illustrated below with a 
sample variable N; 

-0 < 3 - i > 

I, + N. ^ (3.2) 



aN 

a y 


= L-( 

2H V 


N i + 1 - 


a 2 N 




N - 


dy 2 


h A 


. i+1 


a 3 N 


i 


[Y m 


. 3 


„„3 


A i+2 


a y 


2H 



where i is the grid index and H the distance between grid points. 

Centered time differences are used for all quantities except those 
involving friction. The frictional terms are computed at time (t-At). 
In all cases the first step was a forward time step. A special time 
step was introduced every 24 time steps to avoid the separation of 
solutions at even and odd time steps. These special time steps used 
the backwards difference method utilized by Matsuno [1966]. Second 
order equations (2.17), (2.18), (2.19) and (2.20) for time tendencies 
are solved by the exact method of Richtmyer [1957]. A modification of 
the Richtmyer method used by Phillips [1956] was necessary for equation 
(2.22) to satisfy the boundary conditions. Equation (2.21) was solved 
by a direct marching process utilizing the boundary conditions. 

The following boundary conditions were established: 

The quantities A through D = 0 



at ay " 0 and at a y " 0 



all at y = 0 and y = W. 
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Following Phillips [1956] the mean and disturbance vorticity are 
set equal to zero at the boundary in order to make the boundary walls 

smooth. 
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IV. INITIAL CONDITIONS 



A natural initial state would be one of constant temperature and 
no motion. The heating, in time, would build up a baroclinic current. 

No disturbance could develop until the vertical shear of the zonal 
current reached a critical value. In order to save computer time, the 
initial state was chosen to consist of a zonal current with a vertical 
shear near the critical value [see Thompson, 1961]. The initial dis- 
turbance has a wave length and a y-structure such that when superimposed 
on a baroclinic current with no horizontal shear it would have maximum 
growth rate. 

The initial conditions for the model are 



A = 200 sin 7 ^- 



B through D = 0 



f HW 6 



F 



o 



TT 




f HW 6 
o 
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V. RESULTS 



The main experiments were conducted using the physical parameters 
listed in Table I. 



-4 

f = 1 x 10 per second 
o r 

O =2.3 meters squared per second squared per centibar 
squared 

@ = 1.67 x 10 ^ per meter per second 

Ap = 50 centibars 

= 40 meters squared per second 

A = 1 x 10 5 meters squared per second 
H 

a =22.5 degrees 
2 -12 

|j. = 2.49 x 10 per meter squared 

H = 200 kilometers 
At =0.5 hours 
k = 2tt/4000 per kilometer 

Table I - Numerical Values of the Physical Constants 
The variation of experiments 2 and 3 from the main experiment is 
tabulated in Table II. 



- 


W (Km) 


Q (Kj per ton per sec) 


Experiment 1 


8000 


4 x 10 3 (y - ^)/w 


Experiment 2 


8000 


8 x 10 3 (y - j^)/W 


Experiment 3 


4000 


8 x 10‘ 3 (y - £)/W 



Table II - Variation of the Physical Parameters 
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The forecast period for all three experiments was 300 days. Figure 
2 depicts a plot of the disturbance energy for the 300 day forecast 
period. The disturbance energy referred to throughout the discussion 
is defined in Appendix C (Eqs . C-10 and C-ll). A main feature of Fig. 

2 is that the disturbance energy does not reach a steady state although 
it appears to be statistically stable after the first 120 days. 

Figure 3 displays a time section of the mean wind for the period 
120-260 days which is the main time period of interest for this experi- 
ment. The general increase of the disturbance energy begins to level 
off near day 140 which corresponds to the initial appearance of the 
southern jet. This result is similar to observations of Riehl, Yeh 
and LaSeur [1950], The period between day 190 and day 200 is of con- 
siderable interest. The disturbance energy has its greatest fluctuations 
during this period and it corresponds to the termination of the southern 
jet. 

Figures 4-7 represent the amplitudes and phases of Y , Y^ and 
in the form Pcos(kx- 6 ) where P is the amplitude and 6 is the phase. 
Figures 4 and 5 represent day 190, prior to the disturbance peak, and 
Figs. 6 and 7 correspond to day 196, after the peak. The necessary 
condition for barotropic instability is that g - — 5 — change sign some- 

3y 

where in the region [Kuo, 1949]. Calculation at day 190 showed that 
this condition was satisfied in the vicinity of the southern jet. It 
is seen in the lower portion of Fig. 4, between y = 1000 km and y = 

2500 km, that the phase tilts of Y , Y™ and w 0 are opposite to the 

mi 4 

shear of u^. These phase relationships indicate barotropic amplifi- 
cation and therefore a transformation of energy from the mean flow to 
the disturbance. However, the northern portion of Fig. 4 shows the 
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Fig. 3 - Time Cross Section of the Mean Zonal Wind at Level 
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Fig. 4 - Phase Relationship of the Vertical Mean Stream 
Function (Y ) , the Stream Function of the Thermal Wind 
(Yx) and the Vertical Motion (w. = - ^2) f or day 

190. 2 ^ 
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Thousands of kilometers 




Fig. 5 - Amplitude of Y and for day 190. 
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same tilt as the shear of which indicates barotropic damping: this 

implies an energy transformation of disturbance kinetic energy to mean 
kinetic energy. Figure 5 portrays the amplitude of Y and Y m ; clearly 
Y m is larger and both fields have a mid- latitude maximum. A secondary 
peak of amplitude in Y T is beginning to form which is a product of the 
barotropic instability. Figure 6 shows that all parameters (V , Y w 0 ) 
have a phase tilt which indicates barotropic damping over the entire 
range. The amplitude of Y m and V corresponding to Fig. 6 are displayed 
in Fig. 7. Of primary note is the distinct double maximum of amplitude 
which is a result of barotropic instability in the southern region. 
Figure 8 depicts the energy transformations over a period which includes 
the maximum disturbance energy. The first item of note is the nearly 
exact correlation between [ P • P * ] and [P*-K f ] which in an energy balance 
would almost completely offset one another thereby leaving a small 
amount of disturbance potential energy change. It should be noted 
that just before the event of maximum disturbance energy, the barotropic 
damping [K f *K] approached a value of zero. This fact permits the dis- 
turbance kinetic energy to increase at a rapid rate due to the large 
value of [ P f • K 1 ] and to continue to increase until [K f *K] attains a 
value large enough to terminate growth. 

Figure 9 depicts the disturbance energy over the forecast period 
for experiment 2. As stated earlier (Table II), this experiment has 
twice the heating rate of experiment 1. A steady state was not 
attained but the disturbance energy became statistically steady. The 
main differences from experiment 1 are that the disturbance energy 
fluctuates more rapidly and that its magnitude is much higher. 
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A depiction of the disturbance energy over the forecast period for 

experiment 3 is contained in Fig. 10. This experiment differs from 

experiment 1 in that the disturbance between the walls has been cut in 

half. The heating coefficient was doubled from its value in experiment 

1. The disturbance energy reached a steady state condition shortly 

after day 120 and stayed relatively constant to the end of the period. 

Thus we have a single wave traveling through the atmosphere with no 

variation in amplitude. This result is similar to certain dishpan 

experiments with an inner core [Hide, 1953; Fultz, et al, 1959]. 

Figure 11 represents the phases of Y , and w 0 at day 200 and reveals 

m T 2 

that the parameters have similar phase variations. The phase tilts 

indicate barotropic damping since there is a single jet in u^ (see 

Fig. 13). Since Y™ lags Y , there is a northward heat flux and energy 

transformation from the mean potential energy to the disturbance 

potential energy. The correlation mentioned above gives a transformation 

of disturbance potential energy to disturbance kinetic energy since there 
¥ 

is a tendency for cold air to sink and warm air to rise. Figure 12 
represents a plot of the amplitudes of Y^ and y at day 200. The 
amplitude of Y^ is greater than Y^ and both have a single maximum. 

Figure 13 depicts the relationship between u^ , which could be con- 
sidered the surface wind, and u^ at day 200. This mean wind field is 
similar to that observed in the atmosphere, although the gradients are 
much less. At the surface, u is westerly in the mid latitudes and 
easterly at high and low latitudes while at the upper levels the flow 
is westerly with a mid-latitude jet. Figure 14 is a plot of the 
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amplitude of the mean meridional component of the wind (v^) for day 200. 
The amplitude of v, reveals a northerly component at mid latitudes and 
a southerly component at high and low latitudes which implies an indirect 
Ferrel circulation at mid latitudes and a direct Hadley circulation at 
high and low latitudes. 
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Fig. 13 - The Amplitude of the Zonal Component of the mean (u^) 
and the Mean Surface Wind (u^) for day 200. 
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VI. CONCLUSIONS 



A two-level quas i-geos trophic model was formulated for simplified 
studies of the general circulation of the atmosphere. The model was 
further simplified by restricting the disturbance to one wave in the 
zonal direction. 

Numerical integration was computationally stable for a period of at 
least 300 days. In all the experiments the same heating function was 
used which was linear in y and independent of time. The main experi- 
ment with a wall separation of 8000 km produced fluctuating solutions 
with irregular time variations. Attention was focused on a sub period 
of the integration where the disturbance featured a rapid growth followed 
by a rapid decay. Apparently a double jet structure in the mean flow 
became barotropica lly Unstable at the beginning of . the period which 
resulted in the rapid growth of the disturbance. By the end of the 
period, a single jet had replaced the double jet structure. Energy 
transformations verify the roll of [K 1 *K] ,in the evolution of this 
disturbance . 

Experiment 2 featured a double heating rate and results were 
similar to experiment 1 with increased fluctuations and higher energy 
values of the disturbance energy. 

Experiment 3 differed from experiment 1 in that the distance 
between walls was reduced to 4000 km and the heating rate was doubled. 

The disturbance energy became constant after approximately 120 days 
and remained so throughout the forecast period. This situation 
represents a wave of constant amplitude propagating through the 
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atmosphere. The behavior resembles the flow observed in certain dish- 
pan experiments which contained a central core [Hide, 1953; Fultz, et al, 
1959]. Apparently the explanation of the difference in behavior between 
experiment 1 and experiment 3 is related to the fact that for W - 8000 km, 
at least two energy producing modes of disturbance energy can exist in y, 
whereas for W ~ 4000 km only one energy producing mode can exist. The 
presence of two energy producing modes apparently leads to fluctuations 
in time but the one modal case produces a steady wave. 

Further studies utilizing this simple model could be performed such 
as longer period forecasting; changes in the form of heating; simulating 
dishpan experiments by setting 3=0 and concentrating the heating near 
the boundaries, and a change in the wave number k to name a few. Any 
or all of these studies could further enhance our understanding of the 
general circulation of the atmosphere 
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APPENDIX A - COMPUTATION OF VERTICAL MOTION 



The computation of vertical motion (o^) is necessary for the energy 
equation. To derive the computational equation for vertical motion, 
start with 



2 Apa Lbt 



r SY T 
Lc 



dY„ 



^ — I — + u m ^- + V 



m 



dY T H Q 



_2] 

2 J’ 



where 



u m f by 

o y 

, dY 

_ m 

V m f 5x 
o 



5Y 



m 



(A-l) 

(A-2) 
(A-3 ) 



Substituting equations (A-2), (A-3), (2.15) and (2.16) into (A-l), 
collecting terms, and again neglecting terms with wave number 2k or 
higher results in 



2 


|cos kx 


FdC k 


f dE _ 


dF 


“2 " Ap CT 


Ldt ‘ f 
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va y D ' 


‘ dy VJ 



. , rdB k rbE b¥ N1 bF " v 2t 

s ln ta + F (j- C - J- A JJ + ^ + -y. 



*Q, 



f Vdy 

o 



(A-4) 



Both the disturbance and the average vertical motion are given by this 
equation, which is the equation used in the model for the computation 
of vertical motion. 
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APPENDIX B - COMPUTATION OF AVERAGE NORTH -SOUTH WIND 



The computation of the average north-south wind (v^) is necessary 
for the energy calculations. The computation for can be derived 
from Fig. R-l by the following 



dp 



+ V 



“V* - 0 



§r + sr'° 



ou <3v 



-2 + _L.O 

Ap dy 



(B-l) 

(B-2) 

(B-3) 



03 ~~ > 03 <“ 

A v. 



03 



Vf A 



A A V A A 



Fig. B-l 



In finite difference form equation B~3 becomes 
r ) L(- 

1 i Ap V 






VI 



(B-4) 



with s 0,0. 
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APPENDIX C - ENERGY TRANSFORMATION EQUATIONS 



The following definitions are utilized in the energy transformation 
equations of Phillips [1956] 



2 r y, - 2 



Mean Potential Energy P = 



tL 

w 



fM 



dy 



(C-l) 



2 r 2 

Disturbance Potential Energy P' = ^ J(Y'^ - Y'-j) dy (C-2) 



— i r — 2 — 2 

Mean Kinetic Energy K - I (u^ + u., ) dy 



(C-3) 



„dY \^ 2 ^dY lx 2 /ST' 2 



Disturbance Kinetic Energy K 1 - ~ {[C”Sr) + (§t) + C“St) + C^r) 1 dy 

(C-4) 



Phillips [1956] derived the following disturbance energy trans- 
formation equations which provide an excellent means for examining 
the extent to which the model contains the physical processes known 
to be important in the atmosphere. 
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The symbolic notation of the form [ P * P 1 ] for example, signifies a 
transformation of energy from one form - the first in the bracket - to 
the second form, The example [p B F*] represents the transformation of 
mean potential energy into the disturbance potential energy. 

Utilizing equations (2.9), (2.10), (2.15) and (2.16), the components 
of the disturbance energy become 
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Similarly, equations (05) through (09) become 
tp-p'] = - / (B€ - AD) ^ dy 
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